# task: https://docs.google.com/document/d/1snU4dXicuPmFz9XjTF8c7nSw0qGdJAgdJeZVrg-NakI
Высказывание прокурора: “Шанс, что у подсудимого была бы именно такая группа крови, если бы он был невиновен -- всего 1%; значит, с вероятностью 99% он виновен...“.
То есть он оценивал $P(has\_rare\_blood|not\_guilty)$, которую посчитал равной априорной вероятности выбрав произвольного человека, попасть в человека с редкой группой крови. Здесь под rare-подразумевается группой крови как у убийцы (она же редкая).
$$P(has\_rare\_blood|not\_guilty) = P(has\_rare\_blood) = \frac{\sum{people\_with\_rare\_blood\_group}}{\sum{all\_people}}=0.01$$И здесь он уже не прав, потому что: $$P(A) = \sum{P(A|B_i)}$$
То есть: $$P(has\_rare\_blood) = P(has\_rare\_blood|not\_guilty) * P(not\_guilty) + P(has\_rare\_blood|guilty) * P(guilty) = P(has\_rare\_blood|not\_guilty) * P(not\_guilty) + 1 * P(guilty) = 0.01$$
Так как мы знаем, что $P(has\_rare\_blood|guilty)=1$, так как достовено известно, что у убийцы редкая группа крови. И следовательно:
$$P(has\_rare\_blood|not\_guilty) = \frac{P(has\_rare\_blood) - P(guilty)}{P(not\_guilty)}$$И $P(has\_rare\_blood|not\_guilty) = P(has\_rare\_blood)$ лишь при предположении, что $\frac{P(guilty)}{1 - P(not\_guilty)} = 0.01$
При этом в суде должна интересовать другая вероятность, а именно: является ли данный человек убийцей (то есть виновен или нет), считая известным факт того, что у него редкая группа крови. Рассмотрим как меняется вероятность виновности человека, после того, как стало известно, что он является обладателем редкой группы крови, применив теорему Байеса
$$P(guilty|has\_rare\_blood) = \frac{P(has\_rare\_blood|guilty)}{P(has\_rare\_blood)} * P(guilty) = \frac{1}{0.01} * P(guilty)$$Довод адвоката состоял в том, что он предположил, кровь могла равновероятно принадлежать любому человеку с такой редкой группой. И таким образом, он просто посчитал значение как отношение. (Хотя даже в такой постановке кажется странным отбрасывать довод, так как даже если нет никаких других улик, и выбор "виновного" равновероятен среди людей с группой такой крови, то 1e-4 гораздо лучше, чем ~1e-6).
В сущности же нужно рассматривать отношение улики к делу имено как отношение априорной вероятности вины к апостериорной (или относительной): $$\frac{P(guilty|has\_rare\_blood)}{P(guilty)} = \frac{1}{0.01} = 100$$
Некоторые формулы: $$P(seak|positive\_test) = \frac{P(positive\_test|seak) * P(seak)}{P(positive\_test)}$$
Где: $$P(positive\_test|seak) = 1 - p_{false\_negative}$$ $$P(positive\_test|!seak) = p_{false\_positive}$$
$$P(seak) = p_{seak} = 0.01$$$$P(positive\_test) = P(positive\_test|seak) * P(seak) + P(positive\_test|!seak) * P(!seak) = (1 - p_{false\_negative}) * p_{seak} + p_{false\_positive} * (1 - p_{seak}) = 0.01 * (1 - p_{false\_negative}) + 0.99 * p_{false\_positive}$$Подстваляя в формулу: $$P(seak|positive\_test) = 0.01 * \frac{1 - p_{false\_negative}}{0.01 * (1 - p_{false\_negative}) + 0.99 * p_{false\_positive}}$$
Аналогично: $$P(seak|!positive\_test) = \frac{P(!positive\_test|seak) * P(seak)}{P(!positive\_test)}$$
Где: $$P(!positive\_test|seak) = p_{false\_negative}$$ $$P(!positive\_test|!seak) = 1 - p_{false\_positive}$$
$$P(!positive\_test) = P(!positive\_test|seak) * P(seak) + P(!positive\_test|!seak) * P(!seak) = p_{false\_negative} * p_{seak} + (1 - p_{false\_positive}) * (1 - p_{seak}) = 0.01 * p_{false\_negative} + 0.99 * (1 - p_{false\_positive})$$И финально: $$P(seak|!positive\_test) = 0.01 * \frac{p_{false\_negative}}{0.01 * p_{false\_negative} + 0.99 * (1 - p_{false\_positive})}$$
def seak_probability_on_positive_test(p_fp, p_fn):
return 0.01 * (1-p_fn)/ (0.01 * (1-p_fn) + 0.99 * p_fp)
def seak_probability_on_negative_test(p_fp, p_fn):
return 0.01 * p_fn/ (0.01 * p_fn + 0.99 * (1- p_fp))
import numpy as np
import pandas as pd
from plotly import express as px
from plotly.subplots import make_subplots
import plotly.graph_objects as go
from scipy import stats
from sklearn import linear_model
import datetime
import itertools
x = np.linspace(0.95, 1., 100)
y_positive_test_increase_p_fp = np.apply_along_axis(lambda x: seak_probability_on_positive_test(1-x, 0.05), 0, x)
y_positive_test_increase_p_fn = np.apply_along_axis(lambda x: seak_probability_on_positive_test(0.05, 1-x), 0, x)
fig = go.Figure()
fig.add_trace(
go.Line(x=x, y=y_positive_test_increase_p_fp, name='Increasing accuracy of false positive')
)
fig.add_trace(
go.Line(x=x, y=y_positive_test_increase_p_fn, name='Increasing accuracy of false negative')
)
fig.update_layout(
title='Chances of being ill if has positive test depending on accurace of false positive/false negative',
xaxis_title="test accuracy",
yaxis_title="chances of being ill",
)
fig.show()
/Users/wunder9l/opt/miniconda3/envs/made-ml/lib/python3.8/site-packages/plotly/graph_objs/_deprecations.py:378: DeprecationWarning: plotly.graph_objs.Line is deprecated. Please replace it with one of the following more specific types - plotly.graph_objs.scatter.Line - plotly.graph_objs.layout.shape.Line - etc.
x = np.linspace(0.95, 1., 100)
y_negative_test_increase_p_fp = np.apply_along_axis(lambda x: seak_probability_on_negative_test(1-x, 0.05), 0, x)
y_negative_test_increase_p_fn = np.apply_along_axis(lambda x: seak_probability_on_negative_test(0.05, 1-x), 0, x)
fig = go.Figure()
fig.add_trace(
go.Line(x=x, y=y_negative_test_increase_p_fp, name='Increasing accuracy of false positive')
)
fig.add_trace(
go.Line(x=x, y=y_negative_test_increase_p_fn, name='Increasing accuracy of false negative')
)
fig.update_layout(
title='Chances of being ill if has negative test depending on accurace of false positive/false negative',
xaxis_title="test accuracy",
yaxis_title="chances of being ill",
)
fig.show()
С точки зрения пользователя, мне бы хотелось уменьшить число ложно положительных срабатываний: меня бы не устраивала вероятность в 16% после получения положительного теста. С другой стороны, получение отрицательного результата теста дают мне шансы быть здоровым около 99,95%, что кажется вполне приемлимым результатом.
df = pd.read_csv('owid-covid-data.csv', parse_dates=['date'])
df.head()
| iso_code | continent | location | date | total_cases | new_cases | new_cases_smoothed | total_deaths | new_deaths | new_deaths_smoothed | ... | gdp_per_capita | extreme_poverty | cardiovasc_death_rate | diabetes_prevalence | female_smokers | male_smokers | handwashing_facilities | hospital_beds_per_thousand | life_expectancy | human_development_index | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | AFG | Asia | Afghanistan | 2020-02-24 | 1.0 | 1.0 | NaN | NaN | NaN | NaN | ... | 1803.987 | NaN | 597.029 | 9.59 | NaN | NaN | 37.746 | 0.5 | 64.83 | 0.511 |
| 1 | AFG | Asia | Afghanistan | 2020-02-25 | 1.0 | 0.0 | NaN | NaN | NaN | NaN | ... | 1803.987 | NaN | 597.029 | 9.59 | NaN | NaN | 37.746 | 0.5 | 64.83 | 0.511 |
| 2 | AFG | Asia | Afghanistan | 2020-02-26 | 1.0 | 0.0 | NaN | NaN | NaN | NaN | ... | 1803.987 | NaN | 597.029 | 9.59 | NaN | NaN | 37.746 | 0.5 | 64.83 | 0.511 |
| 3 | AFG | Asia | Afghanistan | 2020-02-27 | 1.0 | 0.0 | NaN | NaN | NaN | NaN | ... | 1803.987 | NaN | 597.029 | 9.59 | NaN | NaN | 37.746 | 0.5 | 64.83 | 0.511 |
| 4 | AFG | Asia | Afghanistan | 2020-02-28 | 1.0 | 0.0 | NaN | NaN | NaN | NaN | ... | 1803.987 | NaN | 597.029 | 9.59 | NaN | NaN | 37.746 | 0.5 | 64.83 | 0.511 |
5 rows × 59 columns
start_date = pd.Timestamp(2020, 3, 3)
data_df = df[(df.location == 'Russia') & (df.date >= start_date)]
data_df[['total_cases', 'new_cases']] = data_df[['total_cases', 'new_cases']].astype('int64')
data_df['new_cases'] = data_df['new_cases'].replace(0, 1)
train_df = data_df.head(50)
/Users/wunder9l/opt/miniconda3/envs/made-ml/lib/python3.8/site-packages/pandas/core/frame.py:2963: SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame. Try using .loc[row_indexer,col_indexer] = value instead See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy <ipython-input-9-4438433d6578>:4: SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame. Try using .loc[row_indexer,col_indexer] = value instead See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
pd.options.plotting.backend = "plotly"
train_df.plot(x='date', y=['total_cases', 'new_cases'])
def add_bias(x):
return np.vstack([x, np.ones(len(x))]).T
Используя линейную регрессию, обучите модель с экспоненциальным ростом числа заболевших: y ~ exp(линейная функция от x), где x — номер текущего дня.
actual_x = (train_df.date.values - np.datetime64('2020-03-03'))/ np.timedelta64(1, 'D')
x_for_plot = np.linspace(0, 49, 1000)
X = add_bias(actual_x)
fig = make_subplots(
rows=2, cols=2,
subplot_titles=("Total cases", "Total cases (log)", "New cases", "New cases (log)"),
x_title='Day offset from 03.03.2020',
y_title="Number of cases",
)
y_total = train_df.total_cases.values
model = linear_model.LinearRegression().fit(X, np.log(y_total))
print('Coefs for total_cases: ', model.coef_)
fig.add_trace(
go.Scatter(x=x_for_plot, y=np.exp(model.predict(add_bias(x_for_plot))), name='Prediction curve'),
row=1, col=1
)
fig.add_trace(
go.Scatter(x=actual_x, y=y_total, mode='markers+lines', name='Actual data'),
row=1, col=1
)
fig.add_trace(
go.Scatter(x=x_for_plot, y=model.predict(add_bias(x_for_plot)), name='Prediction curve (log)'),
row=1, col=2
)
fig.add_trace(
go.Scatter(x=actual_x, y=np.log(y_total), mode='markers+lines', name='Actual data (log)'),
row=1, col=2
)
y_new = train_df.new_cases.values
model = linear_model.LinearRegression().fit(X, np.log(y_new))
print('Coefs for new_cases: ', model.coef_)
fig.add_trace(
go.Scatter(x=x_for_plot, y=np.exp(model.predict(add_bias(x_for_plot))), name='Prediction curve'),
row=2, col=1
)
fig.add_trace(
go.Scatter(x=actual_x, y=y_new, mode='markers+lines', name='Actual data'),
row=2, col=1
)
fig.add_trace(
go.Scatter(x=x_for_plot, y=model.predict(add_bias(x_for_plot)), name='Prediction curve (log)'),
row=2, col=2
)
fig.add_trace(
go.Scatter(x=actual_x, y=np.log(y_new), mode='markers+lines', name='Actual data (log)'),
row=2, col=2
)
fig.show()
Coefs for total_cases: [0.19829091 0. ] Coefs for new_cases: [0.18596309 0. ]
Найдите апостериорное распределение параметров этой модели для достаточно широкого априорного распределения. Требующееся для этого значение дисперсии шума в данных оцените, исходя из вашей же максимальной апостериорной модели (это фактически первый шаг эмпирического Байеса).
Вместо обучения целевой переменной $y$ - числа заболевших, будем обучать $z=ln (y)$, которая очевидно связана с целевой перменной. Кроме того, оценка на шум скорее выполняется именно для $z$ (о чем например, свидетельсвтуют графики функций, представленные выше.
В качестве априорного распределения возьмем также нормальное (сопряженно-априорное): $$p(w) = N(\mu_0, \sigma_0)$$ Из графиков выше понятно, что $w_0$ должно быть в районе $[0.2-0.3]$. Однако в задании было сказано выбрать достаточно широкое априорное распределение, поэтому можем взять, например $\mu_0=0, \sigma_0=1$:
ws = [0, 0.5]
sigmas = [0.3, 0.5, 1]
x = np.linspace(-2, 2, 400)
# Create traces
fig = go.Figure(layout=dict(title='Choosing w0, sigma0'))
for w, sigma in itertools.product(ws, sigmas):
fig.add_trace(go.Scatter(x=x, y=stats.norm(w, sigma).pdf(x),
mode='lines',
name=f'w={w}, sigma={sigma}'))
fig.show()
Соответственно, прогнозируемая величина $z_i = ln (y_i)$. Будем исходить из предположения, что каждый y_i независим от предыдущих (хотя очевидно, что это не так). Тем не менее, оценим правдоподобие как: $$LH = p(data|w) = p(\overline{z}|w) = p(z_1,z_2,...z_n|w) = \prod_{i=1}^{n}p(z_i|w)$$ $$ln(p(data|w)) = \sum_{i=1}^{n}ln(p(z_i|w))$$
Будем исходить из предположения, что $z_i ~= f(x_i, w) + noise$ и при этом шум можно распределен нормально относительно предсказанной величины: $$z_i = f(x_i, w) + N(0, \Sigma) = N(f(x_i, w), \Sigma) = N(w^Tx, \Sigma)$$
Тогда апостериорная вероятность также является гауссианой: $$p(w|data) = \frac{p(data|w) p(w)}{p(data)} = C * p(data|w) p(w)$$ $$ln(p(w|data)) = + \sum_{i=1}^{n}ln(p(z_i|w)) = C' -\frac{nd}{2}ln(2\pi) - \frac{n}{2}ln(|\Sigma|) -\frac{1}{2}(X^Tw - \overline{z})^T\Sigma^{-1}(X^Tw - \overline{z}) - ln(2\pi) - \frac12ln(|\Sigma_0|) - \frac12(w-\mu_0)^T\Sigma_0^{-1}(w-\mu_0)=$$ $$= C'' - \frac12w^T(X\Sigma^{-1}X^T + \Sigma_0^{-1})w + (\overline{z}^T\Sigma^{-1}X^T + \mu_0^T\Sigma_0^{-1})w - \frac{n}{2}ln(|\Sigma|) - \frac12\overline{z}^T\Sigma^{-1}\overline{z}$$
Откуда собственно можно получить параметры апостериорного распределения: $$\Sigma_{posterior} = X\Sigma^{-1}X^T + \Sigma_0^{-1}$$ $$\mu_{posterior} = \Sigma_{posterior}(\overline{z}^T\Sigma^{-1}X^T + \mu_0^T\Sigma_0^{-1})$$
$$p(w|data) = N(\mu_{posterior}, \Sigma_{posterior})$$Будем считать, что каждый сэмпл имеет одинаковый и независимый шум. В таком случае: $\Sigma = \sigma^2 I$
И в таком случае легко вычислить соответствующий максимальному правдоподобию шум: $$\frac{d(ln(p(w|data)))}{d\sigma^2} = \frac1{2\sigma^4}(X^Tw - \overline{z})^T(X^Tw - \overline{z}) - \frac{n}{2\sigma^2} = 0$$ $$\sigma^2_{max\_likelihood} = \frac1n(X^Tw - \overline{z})^T(X^Tw - \overline{z})$$
def sigma_max_likelyhood(x, y, w):
err = y - np.matmul(x, w)
return np.mean(err * err)
def bayesian_update(mu, sigma, x, y, sigma_noise):
sigma_inv = np.linalg.inv(sigma)
sigma_n = np.linalg.inv(sigma_inv + (1 / (sigma_noise ** 2)) * np.matmul(np.transpose(x), x) )
mu_n = np.matmul(sigma_n, np.matmul(sigma_inv, mu.T) + (1 / (sigma_noise ** 2)) * np.matmul(x.T, y))
return mu_n, sigma_n
# plot functions
def make_sample_lines(mu, sigma, y_func, n=20, xs=None):
# Посэмплируем и порисуем прямые
my_w = np.random.multivariate_normal(mu, sigma, n)
lines = []
xs = xs or np.linspace(-3, 3, 300)
for w in my_w:
lines.append(go.Scatter(x=xs, y=y_func(w, xs), line=dict(width=0.7),
name=f'w0={w[0]:.2f}, w1={w[1]:.2f}'))
return lines
def make_heatmap(mu, sigma, x=None, y=None):
N = 200
offset = 3
x_min, x_max = mu[0] - max(offset*sigma[0][0], 0.1), mu[0] + max(offset*sigma[0][0], 0.1)
y_min, y_max = mu[1] - max(offset*sigma[1][1], 0.1), mu[1] + max(offset*sigma[1][1], 0.1)
x, y = np.mgrid[x_min:x_max:(x_max-x_min)/N, y_min:y_max:(y_max-y_min)/N]
pos = np.empty(x.shape + (2,))
pos[:, :, 0] = x; pos[:, :, 1] = y
z = stats.multivariate_normal.pdf(pos, mean=mu, cov=sigma)
return go.Heatmap(x=np.linspace(x_min, x_max, N), y=np.linspace(y_min, y_max, N), z=z, showlegend=False, showscale=False,)
iterations_to_plot = [0, 1, 2, 5, 10, 20, 40]
subplot_titles = list(itertools.chain(*[[f'W-distribution after {i} data points', f'After {i} data points'] for i in iterations_to_plot]))
fig = make_subplots(
rows=len(iterations_to_plot), cols=2,
subplot_titles=subplot_titles,
)
X = add_bias((train_df.date.values - np.datetime64('2020-03-03'))/ np.timedelta64(1, 'D'))
y = np.log(train_df.total_cases.values)
cur_mu, cur_sigma = np.array([0, 0]), 2*np.array([[1, 0], [0, 1]])
for i in range(50):
x_step, y_step = X[:i], y[:i]
x_last, y_last = x_step[-1:], y_step[-1:]
if i > 0:
sigma_noise = sigma_max_likelyhood(x_step, y_step, cur_mu)
cur_mu, cur_sigma = bayesian_update(cur_mu, cur_sigma, x_last, y_last, sigma_noise)
if i not in iterations_to_plot:
continue
row = iterations_to_plot.index(i) + 1
fig.add_trace(make_heatmap(cur_mu, cur_sigma), row=row, col=1)
for line in make_sample_lines(cur_mu, cur_sigma, y_func=lambda w,x: w[1]*x + w[0], n=100):
fig.add_trace(line, row=row, col=2)
fig.update_layout(showlegend=False, title_text="Plots on some steps of learning",
height=len(iterations_to_plot)*400, width=400*2)
fig.show()